function [B,G] = biharmonic_bc(dofs,d,g,gn,Bdr,cr,caseNum);
% this is the boundary condition
global V posV T posT E posE TE ET;

m = 0;  % the total number of nonzero entrys
rows = 0;
% find the active boundary:
bdr = find(ET(1:posE,2)==0);
nbdr = length(bdr);
for k=1:nbdr
    tri = ET(bdr(k),1);
    rows = rows + 2*d(tri) + 1;
end
Indx2 = zeros(rows,1);
G = zeros(rows,1);
% must pay attentation to the index turn of A,especially 1,3 and 3,1
A=[2,3;1,3;1,2;3,2;3,1;2,1];
row_begin = 1;
% begin to treat boundary condition one by one bound edge
for k = 1:nbdr
    eg = bdr(k);
    tri = ET(eg,1);
    v1_loc = find(TE(tri,:)==eg);
    v2_loc = mod(v1_loc,3)+1;
    v3_loc = mod(v2_loc,3)+1;
    V1 = V(T(tri,v1_loc),:); 
    V2 = V(T(tri,v2_loc),:);  
    V3 = V(T(tri,v3_loc),:);  
    row_idx = row_begin:row_begin + 2*d(tri);
    tri_dofs = get_tri_dof(dofs,tri);
    if v1_loc==2;
        v1_loc = -2;
    end
    line0 = cr_indices(0,d(tri),v1_loc,cr);
    line1 = cr_indices(1,d(tri),v1_loc,cr);
    [c1,c2] = navstk_first_layer(V2,V3,V1,d(tri),g,gn,caseNum); 
    G(row_idx) = [c1;c2];  
    Indx2(row_idx) = tri_dofs([line0;line1]);
    row_begin = row_begin + 2*d(tri)+1;
end
% [Indx2,I] = unique(Indx2);
% G = G(I); 
N = length(Indx2);
ONE = ones(N,1);
Indx1 = [1:N]';
B = sparse(Indx1,Indx2,ONE,N,max(max(dofs)));

function [c1,c2] = navstk_first_layer(V1,V2,V3,d,g,gn,caseNum);
% This function computes the Bnet coefficients for the navstk eqs over the 1st inward layer 
% of the Line [V1,V2] 
% dirichlet condition
I = [0:d]';
J = d - I;
X = (J*V1(1) + I*V2(1))/d;
Y = (J*V1(2) + I*V2(2))/d;
G0 = feval(g,X,Y,caseNum);
% then the neumann condition
I = [0:d-1]';
J = d - 1 - I;
X = (J*V1(1) + I*V2(1))/(d-1);
Y = (J*V1(2) + I*V2(2))/(d-1);
% comput the first layer which is dirichlet problem
c1 = vdm11(d)\G0;
% comput the normal direction of V1V2
r = sqrt((V2(2)-V1(2))^2+(V2(1)-V1(1))^2);
nx = (V2(2)-V1(2))/r;
ny = (V1(1)-V2(1))/r;
% nx = V3(1)-V1(1);
% ny = V3(2)-V1(2);

Gn = feval(gn,X,Y,caseNum);
% G1 =feval(g1,X,Y,caseNum);  % test for known x and y derivatives
% G2 = feval(g2,X,Y,caseNum);
% Gn = -G2*nx + G1*ny;

[n1,n2,n3] = tcord(V1,V2,V3,[nx,ny]);
% % then the first layer B-ordinates
c2 = (1/d*(vdm11(d-1)\Gn)-n1*c1(1:d)-n2*c1(2:d+1))/n3;

% % the following is right
% c2 = 1/d*(vdm11(d-1)\Gn) + c1(1:d);